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' We investigate the formation processes of a sandpile using numerical simulation. We find 

^ ' a new relation between the fluctuation of the motion of the top and the surface state of 

I ^ ' a sandpile. The top moves frequently as particles are fed one by one every time interval 

^■f^ ' T. The time series of the top location has the power spectrum which obeys a power law, 

■ S{f) ~ and its exponent a depends on T and the system size w. The surface state is 

characterized by two time scales; the lifetime of an avalanche, Tq, and the time required to 

cause an avalanche, Tg. The surface state is fluid-like when Ta ^ Tg, and it is solid-like when 

' Ta -^Tg. Our numerical results show that a is a function of Tg/Ta- 
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1. Introduction 

It is known that the state of a granular system changes between solid-like state and fluid- 
. like state. ^"^^ A sandpile in the formation process is a typical system in which the both states 

■ appear. Its state is solid-like when the feed rate of particles to the sandpile is sufficiently small, 

' and the stress in the sandpile localizes in certain particles. ^"^'^^ Contrastingly, the surface of the 

' sandpile is in fluid-like state when the feed rate is large because avalanches occur frequently, 

, and it is reported that the magnitude distribution of avalanches depends on the grain shape 

and the system size.^^"^^^ In particular, the surface state varies locally and with time, and it 

" I 

^ ' is controlled by the feed rate. Elucidation of the transition of the state is one of interesting 

, problems in granular systems. 

In the previous paper, we studied numerically the formation process of a two- 
dimensional sandpile. We found that the power spectrum of the time series of the top location, 
5'(/), obeys generically a power law, S[f) ~ and that the exponent a depends on the 
feed rate. We defined the left(right) mode as the state which avalanches occur mainly on the 
left(right) slope of the sandpile. In the case where the feed rate is large, if we introduce a 
two- valued function K which takes —1 when the left mode appears and 1 when the right 
mode does, the power spectrum of the time series of K obeys a power law, and its exponent 
is equal to a. 
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In this paper, wc investigate the fluctuation of the top location and avalanches in two- 
dimensional and three-dimensional sandpiles in more detail. Avalanches occur on the surface 
of the sandpile, and the surface state is characterized by two time scales; the lifetime of an 
avalanche, T^, and the time required to cause an avalanche, Tg. We find that the surface state 
which changes between fluid- like state and solid-like state and a depend on Tg/Ta for a two- 
dimensional sandpile. We redefine the left(right) mode as the state that almost avalanches in 
a time interval occur on the left (right) slope and observe the continuation of the left or right 
mode for sufficiently long time. We find that the reason of the continuation is that the memory 
of the mode is stored in the shape of the sandpile. In addition, also for a three-dimensional 
sandpile, the power spectrum of the time series of the top location obeys a power law with 
the exponent which depends on the feed rate. 

This paper is organized as follows. In the next section, we describe the simulation method 
and the setup of the system. In Sec. 3, for the two-dimensional sandpile, it is shown that the 
power spectrum of the top location obeys a power law with the exponent a which depends on T 
and w. In Sec. 4, we consider avalanches in the two-dimensional sandpile. In Sec. 5, the relations 
among Ta,Ts and a are discussed. In Sec. 6, wc present the results for the three-dimensional 
sandpile. Sec. 7 is devoted to discussion and summary. 

2. Discrete Element Method 



We numerically simulate the motion of particles using the discrete element method 
(DEM). Particles are circular in a two-dimensional system or spherical in a three- 
dimensional system with the radii uniformly distributed in the range [0.8d, d]. The force of 
gravity acts on every particle, and elastic force, viscous force and coulomb friction affect each 
pair of particles in contact. Let rrii, Ii and denote the weight, the moment of inertia and 
the radius of the ith particle, respectively. The center of mass, Xj, and the angular velocity 
uji of the ith particle obey the following equations of motion. 




(1) 
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where is the Heaviside function, and iijj and Xij are defined as 



(Xj Xj)/ I Xj Xj 



and 



-^ij — ~t~ I '^i 
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respectively. The normal contact force FnT^ij and the tangential contact force Fj-' are calcu- 
lated as follows. We define Fn as 

Fi^ = F}^e{-F}^), (3) 

where 

= -knXij - TJnnij ■ {±i - ±j) , (4) 

commonly in the two and three dimensional systems. The function Q{—F^ij) means that 
particles arc cohcsionlcss. Parameters kn and r]n represent the spring constant and the viscous 
coefficient in the normal direction. 

We employ different definition of F]-' in the two and three dimensional systems. In the 
two-dimensional system, Ff is defined as in the previous paper, -"^^^ 

Fi^ = ktuihij, (5) 

where tij is the tangential vector, and kt is the spring constant in the tangential direction. 
Displacement is given by the integration of the following equation under the condition 
that the ith and jth particles are in contact, that is when \xj — Xj| < rj + Vi. 

= -((x, - X,-) • t,,- + + T^u:J)Q{^i\Fi^\ - (6) 

where ji is the friction coefficient, and u^ is zero when |xj — Xi| > Vj + ri. In the three- 
dimensional system, F^^ is defined as follows. 

F^ if 
liFn otherwise , 



where 



Fj-' = -fet* - -nt (njj X (xj - Xj) -I- TiUJi -I- r^bjj) x iijj, 

* = ^t; / dt'*(0-tKt'), 

l=\ •^*o 

= (jiUJi rjUj) X njj(t') Xj(t') - Xi(t'), 



Time to is the time when the ith and jth particles begin to contact. Parameter T)t is the 
viscous coefficient in the tangential direction. The tangential vectors, ti and t2, are unit 
vectors perpendicular to n^j. 

There are two physical differences for the tangential forces in the two and three dimen- 
sional systems. One is that viscous term is absent in the two-dimensional system, while it is 
present in the three-dimensional system. Although the latter describes more general cases. 
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Table 


I. Variables and Parameters 




2D 


3D 


rrii 


m (2rj/(i)^ 


m{2ri/df 


h 


mjr?/2 




kn[mg/d] 


1.0 X 10^ 


1.0 X 10^ 


r]n[my/djg\ 


1.0 X 10^ 


1.4 X 10^ 


kt[mg/d] 


2.0 X 10^ 


2.5 X 10^ 


r]t[m^/d/g] 





7.2 X 10 




0.5 


0.2 


w 


20d, 40d, 80d, WOd 


30d 



we emphasize consistency with the previous work.^^) The other is that the tangential shear 
has the maximum value in the two-dimensional system, while there is no limit on the shear 
in the three-dimensional system. In the two-dimensional system, we assume that the shear 
beyond the limit does not affect because particles slip. Although these differences influence 
significantly dense systems with strong shear acting continuously, we infer that the influence 
is noncritical for our sandpile systems because the contact time is not very long. 

Parameter values used in our simulation are listed in Table I, and the physical quantities 
are rescaled to be dimensionless where m represents the weight of a particle with radius d. 
The restitution constants in two and three dimensional systems are respectively about 0.3 and 
0.2 with the values in Table I. 

We make a sandpile on a table which has the origin of coordinates at the center. The 
table in the two-dimensional system is illustrated in Fig. 1 (a). It consists of an alignment of 
particles with diameter d on the x-axis, and its length is w. In the three-dimensional system, 
the table is a fiat circular plate with diameter w on the xy plane as shown in Fig. 1(b). It has 
fixed particles with diameter 0.8d on its fringe. The contact force between particles and the 
plate is calculated in the same manner as that between two particles. 

We carry out simulations using a initial sandpile which is large to cover the table. The 
size of a sandpile is kept virtually constant because particles are eliminated if they fall from 
the table. For a three-dimensional initial sandpile, after we feed sufHciently many particles 
and make an initial sandpile, we fix the particles remained in the sandpile for a long time to 
reduce calculation cost. We use Adams-Bashforth method to calculate the time evolution of 
particles, and the time step is (5t = 1.0 x 10~^. 

We feed particles to the sandpile as follows. Particles are dropped one by one from above 
the origin every time interval T whether avalanches occur or not. The height from which 
particles are dropped, H, is measured from the surface in the two-dimensional sandpile, and 
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'a X 



(b) 




X 
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Fig. 1. The two-dimensional sandpile for w = 80d (a) and (b) the three-dimensional sandpile for 



H is fixed to keep the collision impact given by dropped particle constant. Thus, the height 
from the table changes with time. In the three-dimensional sandpile, we fix the height from 
the table, H', for facilitation of experiments with the same setting. 

3. Fluctuation of the Top Location in the Two-Dimensional Sandpile 

The top location is defined as the center of mass of the highest particle in contact with 
others. This quantity indicates the shape of a sandpile, and the top is moved by avalanches. 

we measure the horizontal position of the top, xtop, in the two-dimensional sandpile and 
calculate the power spectrum of the time series of xtop for given parameter set {w,H,T) 
to characterize the motion of the top. For each parameter set, the power spectrum S{f) is 
obtained through the sample average, and the average is over the power spectra of more than 
ten time series with the length lQ'^\/d/g. In the previous paper, only for w = 80d, we found 
that S{f) obeys a power law, S{f) ~ and that the exponent a increases as T decreases, 
while a is independent of H. 

In this paper, we look into the dependence of a on w. Exponent a is calculated by applying 
the least-squares method for the double logarithmic plot of S{f) in the frequency range 
0.0005 < / < 0.01. Fig.2 shows that a tends to decrease with w and increases drastically as T 
decreases in the range T < lOy^d/^ not only when w = 80d but also in cases where w = 20d 
and 40d. For w = 160d, however, the range in which a changes with T is small. In Sec.5, we 
consider again the dependence of a on u; through the relation between a and avalanches. 



w = 30d (b). 
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Fig. 2. Dependence of a on T and w iov H = 20d 
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Fig. 3. The time-space plot of kinetic energy for {w, H, T) = {80d, 20d, 2^/d/g). 



4. Avalanches in the Two-Dimensional Sandpile 

4-1 Avalanches for small T 

To observe avalanches, we show the time-space plot of kinetic energy for small T in Fig. 3. 
The gray scale represents the value of kinetic energy of particles at position x. The kinetic 
energy is high at x = because particles land there. Avalanches continue for a long time on 
the left slope which is the lower half in Fig.3, and the duration time is sufficiently large in 
comparison with T. 

In the previous paper to measure avalanches on each slope, we calculated kinetic energy 
of particles in the left or right half of the sandpile. The kinetic energy in each part defined as 
the range a;>(iora;<— dis denoted by ki or kr- The magnitude relation between ki and kr 
changes with time. We defined K(t) as 

[ —1, otherwise, 

and we called the state that K = 1{K = —1) the right(left) mode. For small T, the exponents 
of the power spectra of xtop and K are approximately equal in a low frequency range. 
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Fig. 4. The time series of ki and kr in the cases where {w,H,T) = {80d, 20d, l^djg) and 
(80d, 20d, d>Q\/d/g) are plotted. The newly-defined mode is shown with At = 50T. 



4-2 Redefinition of the mode 

To investigate the mode in more detail, we redefine the mode which depend on the state of 
avalanches for a long time scale. Because the current definition determined by instantaneous 
magnitude of avalanche, it is not able to operate when avalanches occur intermittently. 

We redefine the mode using K' defined as follows. 

1, iikr{t)>ki{t), 
K'{t) = \ 0, iikr{t) = ki{t), (9) 

-1, \ikr{t) < ki{t). 

The term, K' = 0, is needed in the case where T is large. We define newly the right(left) 
mode as the state that the time for K' = 1{K' = — 1) amounts to more than 'yAt in the time 
[t — At/2,t + At/2] where 0.5 < 7 < 1, and At is a sufficiently large constant in comparison 
with T. In addition, we introduce the competitive mode defined as the state that the fractions 
of K' = 1 and K' = —1 are comparable. 

In Fig.4, we show the time series of ki and kr and the newly-defined mode with 7 = 0.8 
for T = 2y^d/g and with 7 = 0.6 for T = 80^/d/g. Long-lived modes are observed not only 
for small T but also for large T. 

4-3 Continuation of the mode 

We infer that the mode continues because its memory is stored in either of the shape of 
the sandpile or the motion of particles. To determine which of them is more important factor, 
we carry out examinations as follows. We stop adding particles at a time and restart adding 
after waiting until all particles cease, and we calculate the mode before stop and after restart. 
If the same mode tends to appear before stop and after restart, we consider that its memory is 
stored in the shape because the motion of particles ceases before adding particle is restarted. 
Adding particles is restarted at the time when ki and kr decrease to a constant k*. 
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Table II. In the case of L or R before stop 

After\Before L or R 

Same 58% 

C 40% 

Contrary 2% 

total 139 data 

Table III. In the case of C before stop 

After\Before C 

C 72% 

L or R 28% 

total 211 data 




Fig. 5. The power spectrum of xtop for {'w,H,T) = (SOd, lOOd, ^^djg) in the cases (A), (B) and 
Xf = 0. 



This examination is repeated randomly, and the results show that the memory is stored 
in the shape. Tablell shows the results on the fraction of the mode after restart in the case 
of the left(L) or right (R) mode before stop adding, and Tablelll shows that in the case of 
the competitive(C) mode, for {w,H,T) = {80d,2Qd,2y^d/g) where At = SOT, 7 = 0.8, and 
A;* = 1.0 X 10~^mdg. Although the competitive mode appears more frequently than the left 
or right mode in this criterion, the fraction that the same mode appears in Tablell is higher 
than other modes. In addition, also in Tablelll, the fraction of the same mode is significantly 
high. 

4-4 Relation between the fluctuation of the top location and the position of adding particles 
To clarify relations between the motion of the top and avalanches, we change the position 
of adding particles because we infer that the relations change with the position. We choose 
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randomly the horizontal position at which particles are dropped, x/, in a range and carry out 
experiments for the two cases (A) and (B); (A) is the case where particles are fed uniformly 
above the whole table, — ■u;/2 < Xf < w/2, and (B) is the case where we omit a vicinity of 
the top from the rang, Wg < \xf \ < w/2 where We is a constant. In both cases, we calculate 
the power spectrum of xtop for {w,H,T) = {80d,100d,2y^d/g) and {80d, 100d,80^/d/g) and 
compare with that in the case Xf = 0. 

In the cases of Xf = and (A), the exponents of the power spectra are equal in a low 
frequency range. The thick dashed lines in the Fig.5 and Fig. 6 show the power spectra in the 
case (A). For reference, in the case Xf = 0, we plot the power spectra (thin solid lines) and 
the power functions with the exponent which is the same with that of the power spectra (thin 
dashed lines), respectively. 

In addition, for small T, the exponent of the power spectrum depends on whether particles 
are fed near the top or not. The thick solid lines in Fig.5 and Fig. 6 indicate the power spectra 
in the case (B) where Wg = lOd because — lOd < xtop < lOd in the case Xf = 0. Although the 
exponents in the cases (B) and Xf = are almost equal for T = 80^/d/g, the exponents in 
the case (B) is smaller than that in the case = for T = 2yfdfg. 

For small T, we consider that the motion of the top is different in the cases (B) and x/ = 
because avalanches change with the position of adding particle. Avalanches are frequently 
accelerated by the impact of fed particles for small T, and we infer that the probability of 
the acceleration decreases with the distance between the landing position of fed particles and 
the top because the landing position approaches downstream of avalanches. Therefore, the 
motion of the top changes with the distance. Contrastingly, for large T, in both cases = 
and (B), the probability is low because avalanches induced by a fed particle cease before the 
next particle is fed, hence there is no difference in the motion of the top. 
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Fig. 7. The time series of k for {w, H, T) = {IQOd, 2Qd, h^/dfg) 
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Fig. 8. The time series of ki for {w, H, T) = {20d, 20d, f)-sjd/g) 



5. Fluid-Like State in the Surface of a Sandpile 

In this section, we try to obtain a quantitative relation between the states of the surface 
of a sandpile and the exponent of the power spectrum of xtop, ol. We infer that the state is 
characterized by some time scales for avalanches, and that the motion of the top is related to 
the time scales because the top is moved mainly by avalanches. 

The surface state and a depend on not only T but also w. Actually, a for {w,H,T) = 
(160(i, 20d, 5^/ d/g) is smaller than that for smaller w and the same T as shown in Fig.2. As 
shown in Fig. 7 and Fig. 8, ki for w = 160d is clearly smaller than that for w = 20d, and 
such small kinetic energy is characteristic of the solid-like state, while the state for w = 20d 
is fluid-like. 

5. 1 Time scales for avalanches 

The surface state is related to two time scales for avalanches. One is the time required 
to cause an avalanche, and the other is the lifetime of an avalanche. In the case where the 
former is sufficiently larger than the latter, the state is kept solid-like because the time between 
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Fig. 9. We calculate f{w) from the time series of Ni or Nr with the length 2.0 x for 
{H,T) = (20d,80v^). 



avalanches is long. Contrastingly, the state is fluid-like when these time scales are comparable. 

The lifetime of an avalanche, T^, is independent of w. Lifetime is calculated as the 
average of the duration time in which ki or kr is kept larger than a constant ka- The duration 
time of an avalanche is well-defined when the feed rate is small because each avalanche is 
plainly distinguishable. Therefore, we calculate the time scales for large T. Our numerical 
results with ka = 0.05mdg show that is around 5.0^ d/g for w = 20d, iOd, 80d and IGOd 
when {H,T) = (20^,80^/^7^). 

The time required to cause an avalanche, Tg, depends on T and w. Time Tg is defined as 
the time required to accumulate sufficient amount of particles for causing an avalanche. We 
postulate that Tg is proportional to T, and Tg is represented as follows, 

Ts = Tf{w), (10) 

where f{w) is the typical size of an avalanche and defined as the standard deviation of Ni{t) or 
Nj.{t), where Ni{t) and Nr{t) is respectively the number of particles in the left half and right 
half of a sandpile at time t. The left half is defined as the part in the range — 'u;/2>x>— 1.5d, 
and the right half is the part in the range 1.5d < x < w/2. We find that f{w) increases with 
w as shown in Fig9. 

The results show that the fluid-like state of the surface is kept for a long time if T and w 
are small because Ta and Tg are comparable, and that the state is solid-like if T or u; is large 
because Ta ^ T^. 

The exponent a is related to the surface state and a function of Tg/Ta- We assume that a 
depends on the ratio Tg/Ta and rescale the data in Fig.2 by T* = Ta/f{w). The result is shown 
in Fig.lO. However, to judge whether a depends on only T/T*, more elaborate simulation is 
needed to determine a, Ta and Tg more precisely. 
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Fig. 10. The dependence of a on Ts/Ta when iJ = 20rf and Ta = h^djg 




6. Fluctuation of the Top Location and Avalanches in the Three-Dimensional 
Sandpile 

6. 1 Dependence of fluctuation of the top location on T 

In the three-dimensional sandpile, we measure the top location by the cylindrical coor- 
dinates and calculate the power spectrum of the time series of its azimuthal angle (f) where 
— TT < ^ < TT as in the same manner for the two-dimensional sandpile. The power spectrum 
obeys a power law, S{f) ~ in a low frequency range as shown in Fig. 11. Exponent 
is obtained by a least-square fit of the double logarithmic plot of the power spectrum in the 
frequency range 0.0005 < / < 0.01. Dependence of on T is shown in Fig.l2, which is similar 
to Fig.2 for a in the two-dimensional sandpile, although is larger than a. We consider that 
a reason why < a is because w is small in the three-dimensional systems. 
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Fig. 12. Dependence of on T for {w, H') = {30d, 30d). 



6.2 Relation between the motion of the top and avalanches 

We consider the direction of an avalanche projected on the horizontal plain. The direction 
is represented by the average of particle momentum, (pi,P2)) which is defined by the following 
equation, 



where N denotes the number of particles on the table, and vi^i and Vi^2 are x and y directional 
velocities of the ith particle, respectively. We define the direction of an avalanche as the 
azimuthal angle of the vector (pi,P2), 6, where — vr < <7r. 

To characterize the time series of 9, we show its power spectrum in Fig. 11. The power 
spectrum is proportional to that of ^ in a low frequency range. 

7. Discussion 

In our sandpile system and granular flow in a vertical pipe, there are similar relations be- 
tween the exponent of the power spectrum and the phase space volume of each particle. The 
exponent of the power spectrum of the top location depends on the feed rate in the sandpile, 
and the power spectrum of the density wave in the pipe obeys also a power law^^"^^) with the 
exponent which depends on the inflow rate to the pipe.^*^'^^^ If the power spectrum Sd obeys 
a power law, Sd ~ the exponent /? increases with the volume in the phase space where 
each particle able to move freely. The phase space volume is decreased by restraint conditions 
which are different in the sandpile and flow in the pipe. In the sandpile, the volume in kinetic 
momentum space increases with the feed rate. In this case, avalanches occur frequently, and 
the surface state becomes fluid-like. In the pipe, the volume in kinetic momentum space and 
position space are decreased with the inflow rate because clusters appear in the flow. Devel- 
oping these investigation, for granular systems, it is anticipated that the relations between 
the local state, such as fluid-like or solid-like, and the power spectrum in each systems are 




(11) 



i=l 
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We have investigated relations between the fluctuation of the top location and avalanches 
in formation process of a sandpile using numerical simulations. The top location is moved as 
particles are added one by one every time interval T, and its power spectrum S{f) obeys a 
power law, S{f) ~ /"^, in a long time scale. We found that the exponent a decreases with T 
and the system size w. 

In a two-dimensional sandpile, we defined the right (left) mode as the state that avalanches 
occur mainly on the right(left) slope of the sandpile, and we found that the duration time of 
the left or right mode tends to be long compared to T because the memory of the mode is 
stored in the shape of the sandpile. In a three-dimensional sandpile, the direction of avalanches 
in the horizontal plane changes with time, and the power spectra of the top and the direction 
have the same exponent in a low frequency range for small T. 

The surface state of the sandpile and the exponent a depend on the ratio between the 
lifetime of an avalanche, Ta, and the time required to cause an avalanche, Tg. Our numerical 
results show that Ta is a constant independent of w, while Tg increases with w and T. There- 
fore, the state is kept fiuid-like when T and w are small because ~ T^, and it is solid- like 
when w ov T is large because <C T^. The state relates to the exponent a, and we found 
that a is a function of Tg/Ta- 

Acknowledgment 

I appreciate helpful comments with Hisao Hayakawa, Hiroyuki Tomita, Shinji Takesue, 
Mitsusada Sano and So Kitsunezaki. The numerical calculations were carried out on Altix3700 
BX2 at YITP in Kyoto University. 



14/15 



J. Phys. Soc. Jpn. Full Paper 

References 

R. M. Neddcrman: Statics and Kinematics of Granular Materials (Cambridge, Cambridge, 1992) 
H. M. Jaeger, S. R. Nagel, and R. P. Behringer: Rev. Mod. Phys. 68 (1996) 1259. 
L. P. KadanofT: Rev. Mod. Phys. 71 (1999) 435. 
J. Duran: Sands, Powders, and Grains (Springer, New York, 2000) 
T.Poschel and S. Luding: Granular Gasses (Springer, New York, 2001) 
T. Poschel and N. Brilhantov: Granular Gas Dynamics (Springer, New York, 2003) 
J. P. Wittmer, P. Claudin, M. E. Gates, and J.-P. Bouchaud: Nature 382 (1996) 336. 
L. Vanel, D. Howell, D. Clark, R. P. Behringer, and E. Clement: Phys. Rev. E 60 (1999) R5040. 
J. Ceng, D. Howell, E. Longhi, R. P. Behringer, G. Reydellet, L. Vanel, E. Clement, and S. Luding: 
Phys. Rev. Lett. 87 (2001) 035506. 

J. Ceng, E. Longhi, R. P. Behringer, and D. W. Howell: Phys. Rev. E 64 (2001) 060301. 
V. Prette, K. Ghristensen, A. Malthe-S0renssen, J. Feder, T. J0ssang and P. Meakin: Nature 379 
(1996) 49. 

E. Altshuler, O. Ramos, G. Martinez, L. E. Flores, and G. Noda: Phys. Rev. Lett. 86 (2001) 5490. 
N. Yoshioka: Earth, Planets, and Space 55 (2003) 283. 
G. Urabe: J. Phys. Soc. Jpn. 74 (2005) 2475. 
P. A. Cundall and O. D. L. Strack: Gcotechnique 29 (1979) 47. 
G. Peng and H. J. Herrmann: Phys. Rev. E 49 (1994) R1796. 

G. Peng and H. J. Herrmann: Phys. Rev. E 51 (1995) 1745. 

S. Horikawa, A. Nakahara, T. Nakayama, and M. Matsushita: J. Phys. Soc. Jpn. 64 (1995) 1870. 
S. Horikawa, T. Isoda, T. Nakayama, A. Nakahara, and M. Matsushita: Physica A 233 (1996) 699. 
A. Nakahara and T. Isoda: Phys. Rev. E 55 (1997) 4264. 

O. Moriyama, N. Kuroiwa, M. Matsushita, and H. Hayakawa: Phys. Rev. Lett. 80 (1998) 2833. 
O. Moriyama, N. Kuroiwa, T. Isoda, T. Aral, S. Tateda, Y. Yamazaki, and M. Matsushita: in 
TRAFFIC AND GRANULAR FLOW '01, ed. M. Fukui, Y. Sugiyama, M. Schreckenberg, and D. 
E. Wolf (Springer, New York, 2003) p.437. 

Y. Yamazaki, S. Tateda, A. Awazu, T. Arai, O. Moriyama, and M. Matsushita: J. Phys. Soc. Jpn. 
71 (2002) 2859. 

O. Moriyama, N. Kuroiwa, S. Tateda, T. Arai, A. Awazu, Y. Yamazaki, and M. Matsushita: Prog. 
Theor. Phys. Supp. 150 (2003) 136. 

H. Hayakawa and K. Nakanishi: Prog. Theor. Phys. Supp. 130 (1998) 57. 
H. Hayakawa: Phys. Rev. E 72 (2005) 031102. 



15/15 



